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Abstract. For the stationary one-dimensional nonlinear Schrodinger equation (or 
Gross-Pitaevskii equation) nonlinear resonant transmission through a finite number of 
equidistant identical barriers is studied using a (semi-) analytical approach. In addition 
' to the occurrence of bistable transmission peaks known from nonlinear resonant 

transmission through a single quantum well (respectively a double barrier) complicated 
I (looped) structures are observed in the transmission coefficient which can be identified 

' as the result of symmetry breaking similar to the emergence of self-trapping states in 

^ , , double well potentials. Furthermore it is shown that these results are well reproduced 

by a nonlinear oscillator model based on a small number of resonance cigenfunctions 
of the corresponding linear system. 



> 

^ ■ PACS numbers: 03.65.-w, 03.750.Lm, 42.65.Pc 

OO ' 1. Introduction 

o 

o\ . 

O ■ Transport properties of Bose-Einstein condensates (BECs) in (quasi-) one- dimensional 
waveguides are of considerable current interest, both experimentally and theoretically. 
^ . Especially atom-chip experiments are well-suited to study the influence of interatomic 
interaction on transport properties of BECs in waveguides since various waveguide 
geometries can be realized by different methods [HO]. 

A convenient theoretical approach is based on the one-dimensional Gross-Pitaevskii 
equation (GPE) or nonlinear Schrodinger equation (NLSE) 

ih^{x,t) = l- — —+g\^{x,t)\^ + V{x)j ^{X,t) (1) 

which describes the dynamics in a mean-field approximation at low temperatures 
[11-14]. The nonlinear term g\ip{x,t)\'^ models the interaction between the condesate 
particles. Another important application of the NLSE is the propagation of 
electromagnetic waves in nonlinear media (see, e.g., [15, Ch.8]). The ansatz ilj{x,t) = 
exp{—ifit / h) ipi^x) reduces ([H) to the corresponding time-independent NLSE 

'^^ + 9mx)\' + V{x)-f?j^ix)=0 (2) 
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with the chemical potential (i. 

Within this framework resonant transport through single well (respectively double 
barrier) structures has been studied using analytical and numerical approaches [16-18]. 
It was found that due to the nonlinearity of equation ([2]) the barrier transmission 
coefficient in depencence on the chemical potential /i shows bistable resonance peaks 
which can be related to nonlinear metastable (Siegert) resonance states of the barrier 
potential and described by means of a nonlinear generalization of the Lorentzian profile 
occuring in linear transmission problems [17,19]. These studies correspond to the barrier 
tunneling of coherent monochromatic matter waves with a given chemical potential /i 
that are injected into the waveguide from a BEC reservoir. In the articles cited above 
it was shown that the results obtained from the stationary NLSE ([2]) are in excellent 
agreement with numerical solutions of the time-dependent NLSE 

ihij{x,t) = -—tlj"{x,t) + V{x)tP{x,t)+g\ij{x,t)\^iP{x,t)+foexp{-ifit/h)6{x-xo) (3) 
2m 

where the coupling to a reservoir is modeled by the source term /q exp{—ifit/h)6{x — Xo) 
located at some position x = xq on the left hand side of the barrier (i. e. in the upstream 
region), emitting monochromatic matter waves at chemical potential /i. In contrast to 
the linear case these results cannot be straightforwardly used to predict the scattering 
behaviour of an arbitrary wavepacket since the superposition principle is no longer valid. 

On the other hand double well potentials have been considered in a number of 
theoretical and experimental papers (see e. g. [20-28]). In such systems one observes 
the onset of symmetry breaking and the emergence of new solutions in addition to the 
solutions with linear counterpart for a critical value of the nonlinearity. These results 
strongly motivate a study of related effects in the context of resonant transmission 
through structures consisting of more than one well (or more than two barriers, 
respectively) where one expects the occurrence of both bistability and symmetry 
breaking. The limiting case of resonant transport in infinitely extended periodic 
structures has also been of recent interest (see, e. g. [29-31]). The occurrence of looped 
Bloch bands is one of the major effects of nonlinearity in these systems. Transport 
through a finite number of delta barriers was considered in [32], however, focusing on 
different aspects like the superfluidity of the condensate flow. In this paper we thus 
intend to flU a gap by considering nonlinear resonant tunnelling through a flnite sequence 
of n identical equidistant barriers. For the linear Schrodinger equation transmission 
through such a multi-barrier or truncated periodic potential has been investigated 
in a number of theoretical papers motivated both by experiments with multilayered 
semiconductor heterostructures as well as fundamental issues like their relationship to 
infinitely extended periodic potentials. The particular systems treated in the literature 
include analytically solvable potentials like sequences of rectangular [33, 34] or delta- 
function barriers [34-36]. 

In the following we consider resonant transmission for the stationary one- 
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dimensional NLSE ([2]) with the potential 

V{x) = ^\Y,5{x-jd) (4) 

consisting of n identical delta-function barriers with distance d and strength h'^X/m with 
A > 0. 

This paper is organized as follows. In section [2] we have a brief look at the potential 
dl]) the linear limit {g = 0). In section [3] we introduce a semi-analytical method for 
calculating the transmission coefficient for piecewise constant potentials which is applied 
to the case of double, triple, quadruple and quintuple barrier tunnelling in section HI In 
section [5] these results are compared with the predictions of a nonlinear oscillator model. 
Additional material concerning computational details is presented in an appendix. 

2. Multi barrier transmission in the linear limit 




Figure 1. Transmission coefficients |T„p in dependence on the chemical potential fj, 
for n = 2, 3, 4, 5 and the parameters A = 10, d = 2, g = 0. 



In this section we briefly discuss transmission through the barrier potential (jlj) for 
the linear Schrodinger equation, i. e. equation ([2]) with g = 0. Using the transfer matrix 
technique, Griffiths and Taussig [35] have proven that the transmission coefficient for n 
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delta barriers is given by 



(5) 



with k = y/2mn/h, 

z = cos kd + (X/k) sin kd (6) 

and the Chebyshev polynomials Un{z) of the second kind generated by the recurrence 
relation 

Ur,+liz)=2zUr.iz)-Un-liz) (7) 

starting with Uo{z) = 1 and Ui{z) = 2z. Note that for any given value of z with 
—l<z<l equation ([6]) yields infinitely many solutions k (respectively fi). The 
resonant chemical potentials /xr where |T„(/iR)p = 1 can be determined by solving the 
transcendental equation ([6]) with the roots 

= cos(/7r/r;,), l = l,...,n — l (8) 

of the Chebyshev polynomial f/„,_i(z). Since Un{z) has n zeros, |Tip has no resonances 
and the resonances of |T„p with n > 2 occur in groups of multiplicity n — 1. Because 
of l/n = (i//)/(z/n) with u = 1,2, .. . we see from equation ([8]) that any resonance of 
|T„P is also a resonance of |T;^.„p. This result can be understood in an intuitive way by 
decomposing a series oiv-n identical single barriers into a series of v groups consisting of 
n barriers. The simplest example is given by a quadruple barrier that can be decomposed 
into two double barriers. An incoming plane wave with a chemical potential ^^=2,1=1 
that is in resonance with the double barriers remains an incoming plane wave after 
passing the first double barrier so that it can also pass the second double barrier in the 
same manner. Thus the quadruple barrier is transparent at ^1^=2,1=1 = /in=4,z=2- This 
argument still holds in the case of a finite interaction strength g ^ 0. 

By means of a Taylor expansion of the denominator in equation ([5]) the transmission 
coefficient in the vicinity of a resonance with chemical potential fi = fin,i can be written 
as a Lorentzian (cf. e. g. [33]) 

r^V4 



where 



n,l 



(/i - IJ'n,lf 



A d[/„_i(z) dz 
k dz d/i 



(/i - + TIJA 



2sin2(/7r/ri) 



^(^{n + 2)zUn 



\ dz 


-1 


)dil 





(9) 



(10) 
(11) 

which varies more 



is the full width of the peak at half maximum. The factor sin^(/7r/n 
strongly in dependence on / than the term in the brackets, indicates that within a group 
of resonances the peaks in the middle are broader than the peaks at the sides. 

Figure [1] shows the transmission coefficients for n = 2,3,4,5 for the potential (jlj) 
with X = 10, d = 2 where units with h = m = 1 are used as in all figures and numerical 
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calculations in this paper. Only the first two groups of resonances are shown. It can be 
verified that the positions of the two lowest resonance peaks of IT2P coincide with the 
positions of the second and fifth resonance peak respectively of IT4P as predicted above. 

The groups of resonances and the regions of (almost) zero transmittivity in between 
correspond to the energy bands and band gaps of an infinitely extended delta-comb (or 
Kronig- Penney) potential, respectively (see e. g. [35] and references therein). 



3. Transfer map approach 

Because of the nonlinearity of the GPE the transmission coefficient in the interacting 
case (? 7^ can no longer be obtained by the transfer matrix technique. Instead we 
introduce a method which we call the transfer map approach. To this end we make use 
of an amplitude phase decomposition 



i){x) = ^/S{x)exp{i^{x)) (12) 
which yields the relation 

between the density S{x), the derivative of the phsae $'(x) and the density of the total 
probability current jt. For given values of jt and the transfer map is supposed to map 
the density S{x) and its derivative S'{x) given at some position x on the right hand side 
of the barrier region onto the corresponding quantities S{x) and S'{x) at some position 
X on the left hand side of the barrier region. 

In order to obtain the transfer map of the potential (jlj) we first consider the case 
of a constant potential V{x) = Vq in which an analytical solution for the density S{x) 
is given by (cf. [37,38]) 

S{x) = e + ipdii^igx + 6\p) . (14) 
Using the abbreviation u = gx + 6 the derivative of S{x) and its square are given by 

S'{x) = —2g(ppsn{u\p)cn{u\p)dn{u\p) (15) 
and, by means of the addition theorems of the Jacobian elliptic functions [39], 

S'\x) = [(p _ l)dn2(n|p) + (2 - p)dn\u\p) - dn%u\p)] . (16) 

The parameters in ( fT^ must satisfy 

= -gmif/h^ (17) 

fi-Vo = ^ge + ^g^{2-p) (18) 

mjt + (p - l)9¥^^e - 2(/i - Vo)e^ + 2ge' = . (19) 
Equation (fT8|) can be rewritten as 

^\p-l) = ^'+(3e-^-^l^y. (20) 
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Combining equations and (IT91) we obtain a quadratic equation for tp 



+ 3^ _ '"^ <^ + _^ + - '"^ e = (21) 



with the solutions 




^± = -(;--^l±;/($--^V-2-= + ^^^--^ (22) 



for (yf > (+) and (? > (— ), respectively. Suppose the values of S{x) and S\x) are 
known at some position x = x. The the Jacobian elliptic function dn('u|p) can be 
expressed as 

(hi^{u\p) = {S{x)-e)/^ (23) 
with u = Qx + 5. Thus equation f|T6|) at x = 5 can be written as 

= v\p- l){S{x) -e) + i2-pMSix) - ef - {S{x) - ef .(24) 

Agm 

Using equation ( |T8l) to eliminate <y9 in ( 12^ we finally arrive at the cubic equation 

g \AgmS{x) gS{x) g '7 g 

the real solution of which is e. Now that we know e the value of (f follows from equation 
( !22|) . From equations (ITTl) and (ITS!) we furthermore obtain 

2(/i-K))-3(7£ 



g = ^/^grmp / h and p = 2 . (26) 

Now the only quantity that remains to be computed is the phase 6 of the Jacobi 

elliptic function. This can be obtained by numerically solving equation (l23l) at a; = x. 

However it is more efficient to use the addition theorems 

sn{v\p)cn{u\p)dn{u\p) + sn{u\p)cn{v\p)dn{v\p) 
snygx + d\p) = - oT~r\ \ ) 

1 — psn'^[u\p)sn'^[v\p) 

, . , cn(f |p)cn('u|p) — sn{v\p)dn{v\p)sn{u\p)dn{u\p) . 

cn[gx + o\p) = - (28j 

1 — psn'^[u\p)sn'^[v\p) 

dn{v\p)dn{u\p) - psn{v\p)cn{v\p)sn{u\p)cn{u\p) 

an[gx + o\p) — 2f~\ \ 2^ \ \ ^ ' 

1 — psn'^{u\p)sn'^{v\p) 

instead, where v = g{x—x) and u = gx-\-6. In order to apply these addition theorems the 
values of the Jacobian elliptic functions at x = x are required. The Jacobian function 
dn{u\p) is given by equation fl23l) and the remaining functions can be computed via 
cn{u\p) = cos(amn) and sn{u\p) = sin(am£t) where 



am u = ± arcsin ^ y 1 — dn {u\p) /pj (30) 

is the so-called amplitude of the Jacobian elliptic functions. The sign must be chosen 
such that sgn (sin(am£t) cos(amn)) = sgn {— S' (x) / ip) = sgn(sn(M|p)cn(£t|p)dn(n|p)). 
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Thus for given values of /i — Vq and jt equations ([HD, ([15]), (122]) and (12^ - (15Uj) define a 
map 

W^„vb,i.,x-x : {S{x), S\x)) ^ {S{x), S'{x)) (31) 

whicfi we call the transfer map of the constant potential V{x) = Vq. The matching 
condition for the wavefunction ip{x) and its derivative ip'lx) at the position xq of a 
delta potential with the strength h?\/m is given by ip'{xo—) = ■ip{xo+) and ip'^XQ—) = 
ip'{xo+) — 2Xtp{xo). Straightforward algebra shows that these conditions read 

S{xo-) = S{xo+) and S'{xo-) = S'{xo+) - 4XS{xo) (32) 
in terms of S{x) = \ip{x)\'^ and its derivative. This leads to the map 

Vx : {S, S') ^ {S, S' - AXS) (33) 
for the delta potential with strength h^\/m. Thus we obtain the transfer map 

Mn = {VxU^,j„dT~''^x (34) 

of the potential (jl) with Vq = 0. 

For a stationary scattering state the solution in the region x > {n — l)d, i.e. on 
the right hand side of the potential barriers, is given by a plane wave C exp{ikcx) with 
kc = \/lm{^ii — g\C^')lfi and the total current can be expressed as 

3, = \C\^hkclm. (35) 

The probability density and its derivative at the right hand side of the barriers are 
S[{n — l)(i+) = |Cp and S"((n — l)c/+) = 0. The transfer map (!M]) determines the 
probability density and its derivative at the left hand side of the barriers via 

(5(0_),S'(0-))=A^„(|Cp,0) . (36) 

For the parameters considered in this paper the mean-field interaction potential gS{x) 
outside the potential is negligibly small compared to the chemical potential /i so that 
we can write the wavefunction in the region x < as a superposition Aexplikx) + 
bexp{—ikx) of an incoming and an outgoing plane wave. For the wavefunction ip{0—) 
and its derivative ip'{0— at a; = we thus obtain the condition 

2ikA = ij'{0-) +ikij{0-) , (37) 

with the incoming wave amplitude A and the wavenumber k = y/2mfi/h. In order to 
express this condition in terms of 5'(0— ) and S"(0— ) we multiply it by ip*{0—) arriving 
at 

2ik'^*{0-)A = ^V(O-) + iA;|^(0-)|2 . (38) 
Using S'{0-) = ip*{0-)ij'{0-) + ip*' {0-)ip{0-) and 

jt = -ih{r{0-)tP'{0-) - ^*'(0-)^(0-))/(2m) 

to replace the term ip* {0—)ip' (0—) in equation fl38]) we obtain 

2ik^*{0-)A = S'{0-)/2 + i {kc\C\^ + kS{0-)) (39) 
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where we have used equation to replace jf The absolute square of equation 
provides a convenient condition 

4A;25(0-)2|A|2 = 5'2(o_)/4 + {kc\C\^ + fc^(O-))' (40) 

for the wavefunction in the upstream region. For the parameter range considered in this 
paper the approximation kc ^ k = y/2m^/ h can be made since the effective nonlinearity 
is small outside the barrier region (cf. the discussion above). 

Numerically, for given values of fi and A the squared magnitude | C p of the outgoing 
wave amplitude is obtained by solving the system (136|) . fHOj) . This is achieved by 
combining a bisection method with a finite grid for |Cp. The transmission coefficient is 
then given by \T\^ = jt/jin ~ \C\y\A\\ 

The transfer map (I3T1) of the constant potential V{x) = Vq simplifies considerably 
in the special case S'{x) = in which equation fl25|) reads 



{e-S{x)) 



0. (41) 



gS{x) 

Apart from the trivial solution S{x) = e = const, ip = this equation has the solutions 



= [— ^ V \~9 (^2) 

for (7 > (+) and g < (— ) respectively. For g < equation (ITSl) leads to (cf. [18]) 

V = S{x)-e, (43) 

ioT g > it yields S{x) = e + (p{l — p) and finally, together with equation fl26l) . 
^^2(^^_2^_ 
9 

The phase shift is given by 

6 = -gx {g > 0) or 6 = K{j>) - gx{g < 0) (45) 

(cf. [18]). 

In the following section the transfer map approach is applied to double, triple, 
quadruple and quintuple barrier tunnelling. 



4. Nonlinear mult 1— barrier transmission 



4.1. Double barrier 

Using the transfer map approach described in the previous section we compute the 
transmission coefficient |T2p in dependence on the chemical potential fi for the potential 
dl]) with n = 2 barriers, potential strength A = 10, distance d = 2 (cf. figure [1]) and 
an incoming amplitude A = 0.1 which is displayed in figure [2] for several values of 
the interaction parameter g. We obtain the familiar behaviour for nonlinear single 
well/double barrier tunnelling (see e. g. [16,18,19,40]): For g > 0, the peaks are shifted 
to higher chemical potentials due to the repulsive mean-field term g\ilj{x)\'^ in the GPE. 



Multi-barrier resonant tunneling for the one-dimensional NLSE 



9 




1 1.1 1.2 1.3 0.95 1 1.05 1.1 1.15 



Figure 2. Transmission cocfBcient IT2P in dependence on the chemical potential 
for A = 10, d = 2, A = 0.1. Left panel: .g = (solid blue hne), g = 0.005 (black 
dots), g = 0.1 (red dots). Right panel: .g = (solid blue line), g = —0.05 (black dots), 
g = —0.1 (red dots). 

The wavefunctions of the hnear {g = 0) system corresponding to resonant transmission 
|Tp ^ 1 are more strongly affected by the mean-field term than those corresponding 
to off-resonant transmission because they have a greater total norm dx \ip{x)\'^ inside 
the well. Thus the maximum of a peak experiences a stronger shift than its flanks so 
that the peak bends more and more to the right for increasing nonlinearity g, leading 
to bistability. Analogously a peak bends to the left for an attractive interaction g < 0. 
In other words, in some parameter regions there exist states with the same chemical 
potential but with different average numbers of particles inside the well corresponding 
to different values of the transmission coefficient. The transmission coefficient is thus 
subject to a hysteresis effect as the system has a memory given by the average number 
of particles inside the well. As an example we consider a transmission coefficient for 
repulsive nonlinearity as shown in the left panel of figure [2] for g = +0.1 (red). Let 
us assume that the system is initially prepared in a transmission state corresponding 
to a chemical potential fi ~ 1.1 on the left hand side of the bistable region. If the 
chemical potential fi of the incoming matter wave is slowly increased the values of the 
transmission coefficient follow the upper curve in figure [2] until the end of the bistable 
region is reached. Then the transmittivity "drops down" and follows the only existing 
branch. This behaviour has been explicitly demonstrated in a recent numerical study [41] 
for a double Gaussian barrier using the time-dependent description given in equation 

(ED. 

4-2. Triple barrier 

Figure [3] shows the transmission coefficient \T^\^ in dependence on the chemical potential 
/i for the potential (j4]) with n = 3 barriers, potential strength A = 10, distance L = 2 
(cf. figure [1]) and an incoming amplitude A = 0.1 for various positive values of the 
interaction parameter g. For a moderately repulsive g the resonance peaks bend to 
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Figure 3. Transmission coefficient IT3P in dependence on the chemical potential ^ 
for A = 10, d = 2, A = 0.1. Upper left panel: g = (sohd blue hne), g = 0.017 (black 
dots), g = 0.034 (red dots). Upper right panel: g = 0.036 (black dots), g = 0.039 (red 
dots). Lower left panel: g = 0.05 (black dots), g = 0.1 (red dots). Lower right panel: 
g = 0.25 (black dots), g = 0.5 (red dots). 

the right as known for the double barrier. This effect is shghtly weaker for the second 
peak with larger chemical potential fi than for the first one because the respective kinetic 
energy is higher in comparison with the mean-field interaction energy (cf. the discussion 
in [18]). For g ^ 0.035 a narrow loop which is not connected with the other branches 
of the transmission coefficient emerges close to the second resonance peak. When g is 
further increased the second resonance and the loop approach each other. During the 
process the second resonance peak is slightly deformed until the two structures collide 
and finally unite. A similar behaviour could be observed for the transmission coefficient 
of the finite square well considered in [18]. There, the looped structures originate from 
bound states which have been destabilized and thus turned into resonances due to 
repulsive interaction. Here, the loop is formed by solutions without a linear counterpart 
(so called allochtonous solutions). By comparison with the double barrier case we can 
identify the emergence of unconnected loops as an effect of interaction between the first 
two resonances of the system. 

Figure H] displays the squared magnitudes |'?/'(x)p of the wavefunctions 
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Figure 4. Squared magnitudes |?/'(a;)p of wavefunctions corresponding the peak 
maxima in the transmission coefEcient of the triple barrier with A = 10, d ~ 2^ A = 0.1 
and nonhnearity g = 0.036 (of. the black curve in the upper right panel of figure [3]). 
Solid red: first maximum with fi k, 1.12, |Tp « 1. Dashed dotted black: second 
maximum with /i « 1.22, |Tp « 1. Dashed blue: maximum of the looped structure 
with ^ « 1.23, |r|2 « 0.95. 



corresponding to the peak maxima (respectively the looped structure) in the 
transmission coefficient of the triple barrier with A = 10, = 2, A = 0.1 and nonhnearity 
g = 0.036 (cf. the black curve in the upper right panel of figure[3]). The densities 
corresponding to the two autochtonous states with maximum transmission (solid red and 
dashed dotted black) are symmetric whereas the density corresponding to the maximum 
of the allochtonous loop with |Tp < 1 is asymmetric. Hence this state represents an 
example of symmetry breaking in a nonlinear system similar to, e. g., the self-trapping 
states in double-well potentials (see e. g. [20,22,25,26,28]). 

For even higher values of g the transmission peaks in figure [3] bend more and more 
to the right and the second resonance peak develops into a fork (double peak). Note that 
very narrow structures are not always perfectly resolved because our implementation of 
the transfer map approach uses a finite grid for the outgoing amplitude |Cp (see section 

m>. 

An analogous behaviour can be observed in figure [5] which shows the transmission 
coefficient [Tap for the same potential as in figure [3] for negative values of g. Now the 
curves bend to the left and the looped structure appears in the vicinity of the ffist 
resonance peak. 

4-3. Quadruple barrier 

Now we consider the transmission coefficient IT4P in dependence on the chemical 
potential fj, for the potential (jll) with n = 4 barriers, A = 10, 0? = 2 and A = 0.1. Figure 
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Figure 5. Transmission coefficient IT3P in dependence on the chemical potential ^ 
for A = 10, d = 2, A = 0.1. Upper left panel: g = (solid blue line), g = -0.017 (black 
dots), g = -0.028 (red dots). Upper right panel: g = -0.030 (black dots), g = -0.036 
(red dots). Lower left panehp = —0.05 (black dots), g = —0.09 (red dots). Upper right 
panel: g = —0.15 (black dots), g — —0.3 (red dots). 



[6] displays IT4P for various values of g. As in the case of three barriers for positive 
interaction a narrow structure which is not connected with the other branches of the 
transmission coefficient emerges at the right hand side of the first group of resonances. 
In contrast to the triple barrier case, however, the maximum transmittivity within the 
newly created branch is IT4P 1. For higher values of g the newly created branch unites 
with the third resonance peak. After the unification the transmittivity in the vicinity of 
the third resonance peak no longer reaches full transparency which is another difference 
to the case of three barriers. When g is further increased more and more unconnected 
branches emerge. As for the triple barrier the situation is completely analogous in the 
case of negative interaction. By comparing the respective transmission coefficients for 
= and g = ±0.1 in figures [6] and [2] we can directly verify a property predicted in 
section[2], namely that the position of the second resonance peak of the quadruple barrier 
always coincides with the position of the first resonance peak of the double barrier. 
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Figure 6. Transmission coefficient IT4P in dependence on the chemical potential /i 
for A = 10, d = 2, A = 0.1. Upper left panel: .g (sohd blue hne), g = 0.015 (black 
dots), g = 0.02 (red dots). Upper right panel: g — 0.04 (black dots), g = 0.1 (red 
dots). Lower left panel: 5 = (sohd blue hne), g = —0.01 (black dots), g = —0.02 
(red dots) . Lower right panehg — —0.04 (black dots), g — —0.1 (red dots) 



4.4- Quintuple barrier 

Adding another delta function we arrive at the quintuple (n = 5) barrier with the 
parameters A = 10, d = 2 and A = 0.1 (cf. figure [1]). The respective transmission 
coefficient IT5P in dependence on /x is shown in figure [7] for several values of the 
interaction constant g. Similar to the case of the quadruple barrier (figure [6]) increasing 
g leads to the formation of an unconnected structure on the right hand side of the first 
group of resonances, its unification with the highest resonance within this group as well 
as to the emergence of more unconnected branches. In addition, the second resonance 
peak of the group develops into a fork of tree subpeaks. Again, the system reveals a 
completely analogous behaviour for attractive interactions g < 0. 

5. Nonlinear oscillator model 

In [19,40] it was shown that nonlinear resonant tunnelling can be understood in terms of 
Siegert resonances, i.e. in the vicinity of a resonance the wavefunction tlj[x) of the system 




Figure 7. Transmission coefficient IT5P in dependence on the chemical potential fi 
for A = 10, d = 2, A = 0.1. Upper left panel: g^O (sohd blue hne), g = 0.008 (black 
dots), g = 0.013 (red dots). Upper right panel: g = 0.02 (black dots), g = 0.05 (red 
dots). Lower left panel: g = (soHd blue hne), g = —0.02 (black dots), g = —0.01 
(red dots). Lower right panel: g — —0.04 (black dots), g — —0.1 (red dots). 



is approximately given by a so-called skeleton wavefunction ipsk{x) which satifies purely 
outgoing (Siegert) boundary conditions and yields a complex eigenvalue fisk — irsk/2. 
Thus we only take into account the resonant contribution to transmission neglecting the 
non-resonant contribution originating from sequential single barrier tunnelling. 

This approximation can be incorporated in the time-dependent description of 
nonlinear resonant tunnelling mebtioned in the introduction where a source term is used 
to model the injection of an incoming coherent matter wave with chemical potential fi. 
Inserting the ansatz ip{x,t) = exp{—ifit / h)ipsk{x) into equation ([3]) yields 

{Hq- ^ + g\^sk{x)\'^)i'sk{x) + i/o^(a; - xq) = (46) 

where we have chosen a constant source strength f{t) = fo located at some position 
Xo < 0. In analogy to [28] we expand the skeleton wavefunction ipsk{x) in a Galerkin-type 
ansatz 

riB 

ipskix) = ^CjUj{x) (47) 



Multi-barrier resonant tunneling for the one-dimensional NLSE 



15 




1 

0.8 
0.6 
0.4 
0.2 




0.95 1 



o v 



1.05 1.1 1.15 



Figure 8. Transmission coefficient IT2P in dependence on the chemical potential /i for 
A = 10, d = 2, = 0.1. Left panel: g = 0.05. Right panel: g = -0.05. The stability 
predictions of the transfer map approach are indicated by black dots, the results of the 
nonlinear oscillator model by blue asterisks (stable regions) and red circles (unstable 
regions). 



using the first u-q eigenf unctions {uj} and respective eigenvalues {yUj 
linear [g = 0) system 



iTj/2} of the 



(4^ 



with Siegert boundary conditions. 



The eigenfunctions are made square-integrable 

To calculate the 



by means of exterior complex scaling (see appendix Appendix A ), 
stationary states we insert the ansatz Wl\ into equation fj46|) and consider its projections 

2 



Cj(/ij -iTj/2 - ^) + g / dxf 



i=l 



CiUi{x) 



^QMi(x) + ifoV*{xo) = 



(49) 



on the TT-B left eigenvectors {vj} of Hq. The nonlinear equations fH9|l . which determine 
the riB coefficients {cj}, are solved with a Newton algorithm. Obviously all equations 
decouple in the noninteracting case g = 0. A system of nonlinear coupled oscillators 
similar to the one described by equation (1491) has been investigated in the context of 
micromechanical and nanomechanical resonator arrays [42]. 

Since the transmission coefficient for a potential with n barriers shows groups of 
n—1 resonances we use ub = n—1 basis functions to compute the transmission coefficient 
in the vicinity of the first group of resonances. 

The source strength /o is connected with the incoming wave amplitude A via 
/o = h'^ikA/m (cf. [16,19,40]) with k = \J2m^/h. For simplicity we choose Xq = 0. The 
transmission coefficient is given by the solutions of (H9|) via 

|T|^ = — 

Jin 



(50) 



where 



= {n-l)d 



(51) 
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Figure 9. Upper panels: Transmission coefficient jTap in dependence on the chemical 
potential ^ for A = 10, d = 2, A = O.f . Left: g = 0.0366. Right: g = O.l. The results of 
the transfer map approach are indicated by black dots, the stability predictions of the 
nonlinear oscillator model by blue asterisks (stable regions) and red circles (unstable 
regions). Lower panels: corresponding occupation numbers of the ground mode |cip 
(solid blue line) and the first excited mode |c2p (dashed red line). 



and Jin = hk\A\'^/m. In this resonance ansatz the system is described by a small number 
of square integrable functions rather than by a continuum of distributions which can be 
favourable in many situations. Another advantage lies in the fact that the stability of a 
stationary solution can be analyzed in a straightforward way (see below). 

For illustration we have a closer look at the special case of a single mode, i. e. tt-b = 1 
(and thus n = ne + l = 2), which models tunnelling through a single well/double barrier 
structure. Equation ( l49l) now reads 

ci(/ii - iri/2 - fi) + gwll\ci\^ci + ifovKxo) = (52) 

with wll = f'^^vl{x)ul{x)ui{x)ui{x) dx. The squared magnitude of equation (152|) 



l/o 


2 


vi{xo)\'^ 


(/il + gRe{wl\)\ci 


2 _ 


^)2 + (ri/2 + (?Im«)|ci 


2^2 



provides a self-consistent equation for the occupation number |cip of the basis function 
Ui{x). Note that, due to symmetry, |fi(xo)P = |Mi(a;o)P = \ui{d+ \xo\)\'^. By means of 
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Figure 10. Upper panels: Transmission coefRcient jTap in dependence on the chemical 
potential ^ for A = 10, = 2, A = 0.1. Left: g = -0.03. Right: g ^ -0.09. The results 
of the transfer map approach are indicated by black dots, the stability predictions of the 
nonlinear oscillator model by blue asterisks (stable regions) and red circles (unstable 
regions). Lower panels: corresponding occupation numbers of the ground mode |cip 
(solid blue line) and the first excited mode |c2p (dashed red line). 



the Siegert formula |?ii(a;o)| can be expressed in terms of the decay coefficient Vi via 
ri/2 = ,J , ^ ^ xo 54 

with J^^^'^''' dx |mi(x)P ~ 1 and k = y/2mfi/h. In order to express equation (l53l) in 
terms of |Tp = jt/jm instead of |cip we evaluate the current density jt a.t x = d + \xo\ 
which yields jt = ^\ui{xo)\'^\ci\^. Using equation ( 154|) we obtain jt = ^|cipri/2. The 
transmission coefficient is thus given by 

Using equations ([SID, dSS]) and |/oP = h'^k^\A\'^ /m equation (1551) can be written as 

' ' (/. - /..kP + ry4 



(56) 
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with the skeleton curves 

/Xsk(lTp) =/ii+^7ReKl)^^^|rr (57) 

ml 1 

r3k(|Tn/2 = ri/2 + (^ImK|) ^'^''^'V r • (58) 

ml 1 

For small values of Im(wJ| we can make the aproximation Fi Fg^ in the numerator of 
equation (l56i) arriving at the nonlinear Lorentz profile 

r2 I A 

' ' (/^-/^sk)^ + ry4 ^''^ 

derived in [19,40]. Here the skeleton curves fl57j) and fl58l) are given in first order 
approximation in |Tp. 

In order to perform a linear stability analysis of a stationary solution ipsk we insert 
il){t) = {ipsk + ^"^{t)) exp{—ifit / h) into equation ([3]) and retain only terms linear in 
The spectral decomposition 6ip{t) = X-6xp(— icjt) + x^exp(iu;*t) then leads to the 
Bogoliubov-de Gennes equations 

with -f^GP = Hq + g\ipsk\'^- Instability occurs if there are eigenmodes with positive 
imaginary part since their population grows exponentially in time. The eigenvalue 
equation (jHOD is solved in the dual basis {uj} and {vj}. 

In figures [814TT] we compare the predictions of the nonlinear oscillator model with the 
results of the transfer map approach for different numbers of barriers and interaction 
constants. In all cases the agreement between both methods is quite good. Within 
the nonlinear oscillator approach branches which are not connected to the main part 
of the transmission coefficient prove difficult to find numerically and are therefore not 
taken into account. The stability predictions for the double barrier (figure [8]) agree with 
the results expected for a single parametrically driven nonlinear classical oscillator (see 
e.g. [43]) or quantum oscillator in the classical (mean- field) limit [44]. This is also in 
agreement with recent numerical results for a double Gaussian barrier [41] where the 
dynamical stability of the upper branch of the transmission is explicitly demonstrated 
by means of a time-dependent simulation. For more than two barriers the model 
predicts an increasingly complicated distribution of stable and unstable regions. The 
occupation numbers {|cjp} of the modes {uj} shown in the lower panels of figures [9HTT] 
indicate that the autochtonous branches of the respective transmission coefficients are 
mainly described by one mode only whereas the allochtonous branches are formed by 
superpositions of two or more modes. 

6. Conclusion 

In this paper we considered nonlinear resonant tunnelling through sequences of n 
identical and equally spaced delta barriers. The stationary transmission states were 
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Figure 11. Upper panels: Transmission coefficients in dependence on the clremical 
potential for A = 10, = 2, A = 0.1. Left: |T4|2 for g = 0.03. Riglit: [Tsp for 
g = 0.05. The results of the transfer map approach are indicated by black dots, the 
stability predictions of the nonlinear oscillator model by blue asterisks (stable regions) 
and red circles (unstable regions). Lower panels: corresponding occupation numbers 
of the ground mode |cip (thin blue line), first |c2p (dashed red line), second |c3p 
(dashed dotted black line) and third excited mode |c4p (bold green line). 



calculated by means of a transfer mapping approach based on the complex solutions 
of the free time-independent NLSE given by Jacobi elliptic functions. As observed 
for single well/double barrier tunnelling (see [16, 18]) the nonlinearity renders the 
transmission coefficient bistable in the vicinity of a resonance. In addition, looped 
structures appear, which are not connected with other branches of the transmission 
coefficient. If the interaction parameter g is further increased these structures unite 
with the main part of the transmission coefficient through an inverse beak-to-beak 
bifurcation. A similar effect was observed in the transmission coefficient of the finite 
square well (see [18]) for branches of the transmission coefficient originating from bound 
states of the linear {g = 0) system destabilized by interaction. Increasing the number of 
barriers and the nonlinearity leads to the emergence of more and more complicated 
structures in the transmission coefficient which result in a suppression of resonant 
transport. 

Comparison with a finite basis calculation based on the resonance wavefunctions of 
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the linear system shows that the effects described above can be understood in terms of 
nonhnear parametrically driven coupled oscillators. The finite basis approach also offers 
a straightforward way to analyze the stability of different branches of the transmission 
coefficient by solving the corresponding Bogoliubov-de Gennes equations. 
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Appendix A. Left and right resonance eigenfunctions in the linear limit 

The nonlinear oscillator approach in section [5] requires the computation of the resonance 
eigenfunctions u{x) and corresponding eigenvalues /i — ir/2 of the Hamiltonian Hq given 
in equation ( HHl) with the potential V{x) = (^^/?ti)A ^"~q — jd) given in equation 
(jl]) which are obtained by solving the stationary Schrodinger equation 

-^dl + V{x) ^ u{x) = (/i - \T/2)u{x) (A.l) 

with Siegert boundary conditions. We make the ansatz 

{exp(— i/cx) X < 
IjSin{kjd + ^j) {j -l)d<x <jd,0 <j <n-l (A. 2) 

exp(ifcx) X > {n — l)d 

with k = ^2m(/i — iT/2)/h which satisfies the Siegert boundary conditions 
lim^^^-i-oo u'{x) = ±ik u{x). The matching conditions at x = and x = {n — l)d read 

1 = Jl sin(i9i) , -ik = kh cos(i9i) - 2A (A.3) 

and 

/c/„_i cos(?9„_i) =ik-2\. (A.4) 
At X = jd, < j < n — 1 we obtain 

Ij sm{kjd + 'dj) = Ij+i sm{kjd + ^j+i) , (A. 5) 

klj cos{kjd + i^j) = klj+i cos{kjd + — sm{kjd + ^j) ■ (A. 6) 

These equations are solved numerically for the complex quantities k, 'dj and Ij, 

< j <n. 

The wave function u{x) diverges for x — > oo since Im(fc) < 0. Therefore we use 
exterior complex scaling (see e. g. [45,46]) to make the wave function square integrable. 
The X coordinate is rotated by an angle 6^ from the point where the potential V{x) 
becomes zero. In our case this reads 

xexp(i^c) X < 

X ^ { X X < < (n - . (A.7) 

(n — l)d + (x — (n — l)d) exp(i6'c) x > {n — l)d 
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In the scaled region the Schrodinger equation becomes exp(2i6'c)M"(x) + = 0. The 

matching conditions flA.3|) and flA.4l) remain unaltered. For a sufficiently large rotation 
angle 9c the wavefunction u{x) becomes square integrable in < x < oo. 

Since Hq is symmetric, the corresponding left eigenfunctions v{x) are given by 
= We normalize the eigenstates such that (ixv*{x)u{x) = 1. 
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